// -*- mode: c++; indent-tabs-mode: nil; -*-
//
//The Biomolecule Toolkit (BTK) is a C++ library for use in the
//modeling, analysis, and design of biological macromolecules.
//Copyright (C) 2001-2006, Eric Alm <ealm3141@users.sourceforge.net>,
//                         Tim Robertson <kid50@users.sourceforge.net>
//
//This program is free software; you can redistribute it and/or modify
//it under the terms of the GNU Lesser General Public License as published
//by the Free Software Foundation; either version 2.1 of the License, or (at
//your option) any later version.
//
//This program is distributed in the hope that it will be useful,  but
//WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//Lesser General Public License for more details.
//
//You should have received a copy of the GNU Lesser General Public License
//along with this program; if not, write to the Free Software Foundation,
//Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

// contact_order.c

#include <iostream>
#include <cmath>
#include <cstdlib>

#include <btk/core/math/linear_algebra.hpp>
#include <btk/core/atoms/pdb_atom.hpp>
#include <btk/core/molecules/monomer.hpp>
#include <btk/core/molecules/polymer.hpp>
#include <btk/core/io/pdb_system.hpp>

typedef BTK::ATOMS::PDBAtom atom;
typedef BTK::MOLECULES::Monomer<atom> monomer;
typedef BTK::MOLECULES::Polymer<monomer> polymer;
typedef BTK::IO::PDBSystem<polymer> pdb_system;

using std::cout;
using std::endl;
using BTK::MATH::within_sqr_dist;

int main(int argc, char * argv[])
{
  double cutoff = 6.0;
  std::string filename;

  if (argc == 2) {
    filename = argv[1];
  } else if (argc == 3) {
    filename = argv[1];
    cutoff = std::atof(argv[2]);
  } else {
    std::cerr << "usage: " << argv[0] 
              << " <pdb filename> [distance cutoff (default = 6.0)]" 
              << std::endl;
    return -1;
  }

  pdb_system pdb(filename);

  pdb_system::atom_iterator ai, aj;
  int counts = 0;
  double order = 0;
  double cutoff_sqrd = cutoff*cutoff;

  for (ai = pdb.structure_begin(); ai != pdb.structure_end(); ++ai) {
    for (aj = pdb.structure_begin(); aj != pdb.structure_end(); ++aj) {
      if (ai->res_number() == aj->res_number()) continue;
      
      if (within_sqr_dist(ai->position(),aj->position(),cutoff_sqrd)) {
        order += std::abs(ai->res_number() - aj->res_number());
        ++counts;
      }
    }
  }

  double contact_order = order/counts;
  
  cout << "Contact Order: " << contact_order << endl;
  cout << "Relative Contact Order: " << contact_order/pdb.num_monomers() << endl;

  return 0;
}



